#################################################################
##DO SELF-REPORTING REGIMES MATTER? EVIDENCE FROM THE CAT########
## Cosette D. Creamer and Beth A. Simmons #######################
## February 2019 ###############################################
### Replication Script for IPTW MSM Analysis with CIRI ##########
###############and CIRI-LHRPS Interaction (Appendix B)###########
#################################################################
#################################################################

## Clear environment
rm(list=ls())
sink("CAT CIRI MSM Analysis.txt")

### Libraries
library(doParallel)
library(foreach)
library(stargazer)

number_core<-detectCores()-1
cl<-makeCluster(number_core)
registerDoParallel(cl)

##################################################

foreach (alphavalue =seq(0,0, 0.01)) %dopar% {
  iterations <-20000
  mod<-seq(1,3, 1)
  for (m in mod){
    model <- m
    lead <- 2
    seed <- 1234
    alpha<-alphavalue  
    if (!is.na(seed)){
      set.seed(seed, kind="L'Ecuyer-CMRG")
    }
    
    
    ##################################################
    ####### Setting Up R     					  		  ########
    ##################################################
    
    ### Libraries
    library(data.table)
    library(MASS)
    library(car)
    library(SparseM)
    library(ggplot2)
    library(survey)
   
    
    ##############################################
    ### Data loading and pre-processing
    ##############################################
    
   
    ## Load CAT data (full)
    sink("CAT CIRI MSM Analysis.txt",append=TRUE)
    catfull <- read.csv("CATReporting_Final.csv", header=TRUE)
    
    ## Subset data to only include observations for CAT ratifiers, for observations > t+1 (t = year of ratification)
    cat <- subset(catfull, catfull$CATob==1)
    
    ### If we're running model 1
    if (model == 1){
      ## Remove observations with missingness for covariates used - REPORTING
      
      cat <- subset(cat, !is.na(cat$gdpl1) & !is.na(cat$poptot) & !is.na(cat$tortl1) & 
                      !is.na(cat$polityl1) & !is.na(cat$CATreview)  & !is.na(cat$nhril1) 
                    & !is.na(cat$CATart22l1)  & !is.na(cat$CATregiondensityl1) & !is.na(cat$nevdem) & 
                      !is.na(cat$regionall1) & !is.na(cat$transdeml1) & !is.na(cat$hrtprop) & !is.na(cat$transition2l1))
      
      
      cat$CATregiondensityl1 <- 100*cat$CATregiondensityl1
      cat$HRTpercentl1 <- 100*cat$hrtpropl1
      cat$yr <-as.factor(cat$year)
      cat$lngdpl1 <- log(cat$gdpl1)
      cat$lnpop <- log(cat$poptot)
      
      cat$tortf1<-factor(cat$tortf1, ordered=TRUE)
      cat$tortf2<-factor(cat$tortf2, ordered=TRUE)
      
      
      ## Define a dataset with no missing obs for torture t+1 or t+2 (depending on argument)
      if (lead == 1){
        catuse <- subset(cat, !is.na(cat$tortf1))
        catuse$outcome <- catuse$tortf1
        catuse$outcome<-factor(catuse$outcome, ordered=TRUE)
      }else if (lead == 2){
        catuse <- subset(cat, !is.na(cat$tortf2))
        catuse$outcome <- catuse$tortf2
        catuse$outcome<-factor(catuse$outcome, ordered=TRUE)
      }else{
        stop("'lead' argument must be 1 or 2")
      }
      
      ##### Sort dataset by country- year
      catuse <- catuse[order(catuse$cowcode, catuse$year),]
    }else if (model == 2){
      cat <- subset(cat, !is.na(cat$gdpl1) & !is.na(cat$poptot) & !is.na(cat$tortl1) & 
                      !is.na(cat$polityl1) & !is.na(cat$CATreview)  & !is.na(cat$nhril1) 
                    & !is.na(cat$CATart22l1)  & !is.na(cat$CATregiondensityl1) & !is.na(cat$nevdem) &
                      !is.na(cat$regionall1) & !is.na(cat$transdeml1) & !is.na(cat$transition2l1) 
                    &!is.na(cat$hrtpropl1) & !is.na(cat$transparencyindex_l1))
      
      cat$CATregiondensityl1 <- 100*cat$CATregiondensityl1
      cat$HRTpercentl1 <- 100*cat$hrtpropl1
      cat$yr <-as.factor(cat$year)
      cat$lngdpl1 <- log(cat$gdpl1)
      cat$lnpop <- log(cat$poptot)
      
      cat$tortf1<-factor(cat$tortf1, ordered=TRUE)
      cat$tortf2<-factor(cat$tortf2, ordered=TRUE)
      
      
      ## Define a dataset with no missing obs for torture t+1 or t+2 (depending on argument)
      if (lead == 1){
        catuse <- subset(cat, !is.na(cat$tortf1))
        catuse$outcome <- catuse$tortf1
        catuse$outcome<-factor(catuse$outcome, ordered=TRUE)
      }else if (lead == 2){
        catuse <- subset(cat, !is.na(cat$tortf2))
        catuse$outcome <- catuse$tortf2
        catuse$outcome<-factor(catuse$outcome, ordered=TRUE)
      }else{
        stop("'lead' argument must be 1 or 2")
      }
      
      ##### Sort dataset by country- year
      catuse <- catuse[order(catuse$cowcode, catuse$year),]
      ### Interaction model
    }else if (model == 3){
      
      cat <- subset(cat, !is.na(cat$gdpl1) & !is.na(cat$poptot) & !is.na(cat$tortl1) & !is.na(cat$farissl2) & 
                      !is.na(cat$polityl1) & !is.na(cat$CATreview)  & !is.na(cat$nhril1) 
                    & !is.na(cat$CATart22l1)  & !is.na(cat$CATregiondensityl1) & !is.na(cat$nevdem) &
                      !is.na(cat$regionall1) & !is.na(cat$transdeml1) & 
                      !is.na(cat$transition2l1) & !is.na(cat$transparencyindex_l1) &!is.na(cat$hrtpropl1))
      
      cat$CATregiondensityl1 <- 100*cat$CATregiondensityl1
      cat$HRTpercentl1 <- 100*cat$hrtpropl1
      cat$yr <-as.factor(cat$year)
      cat$lngdpl1 <- log(cat$gdpl1)
      cat$lnpop <- log(cat$poptot)
      
      cat$tortf1<-factor(cat$tortf1, ordered=TRUE)
      cat$tortf2<-factor(cat$tortf2, ordered=TRUE)
      
      
      ## Define a dataset with no missing obs for torture t+1 or t+2 (depending on argument)
      if (lead == 1){
        catuse <- subset(cat, !is.na(cat$tortf1))
        catuse$outcome <- catuse$tortf1
        catuse$outcome<-factor(catuse$outcome, ordered=TRUE)
      }else if (lead == 2){
        catuse <- subset(cat, !is.na(cat$tortf2))
        catuse$outcome <- catuse$tortf2
        catuse$outcome<-factor(catuse$outcome, ordered=TRUE)
      }else{
        stop("'lead' argument must be 1 or 2")
      }
      
      ##### Sort dataset by country- year
      catuse <- catuse[order(catuse$cowcode, catuse$year),]
      
    }else{
      stop("Model must be 1,2, or 3")
    }
    
    ###################################################
    ##### Main estimation/bootstrap routine for model 1
    ###################################################
    
    ## Reporting model (alone)
    if (model == 1){
      ## Calculate number of country clusters
      n_country <- length(unique(catuse$cowcode))
      
      ## Number of coefficients
      n_coef <- 5
      
      ## Placeholder for holding bootstrap coefficients
      boot_coefs <- as.data.frame(matrix(nrow=iterations, ncol=(n_coef+2)))
      
      # Start bootstrap loop - since we're concerned about sampling datasets that have non-convergence
      # we make this a while loop and use tryCatch when dealing with non-feasible datasets.
      
      # Initialize iteration number
      iter <- 1
      while(iter <= iterations){
        print(iter)
        # If it's not the first run, bootstrap
        if (iter != 1){
          # Resample dataset by country clusters
          boot.data <- data.frame()
          for (kit in 1:n_country){
            cnum <- sample(unique(catuse$cowcode), size=1) # Sample a country
            sampled.country <- catuse[catuse$cowcode == cnum,]
            sampled.country$uniqueID <- kit
            boot.data <- rbind(boot.data, sampled.country) # Append to the bootstrap dataset
          }
          # If it's the first run, use the raw data
        }else{
          boot.data <- catuse
          boot.data$uniqueID <- boot.data$cowcode
        }
        
        ### Get Coefficients - catch any errors
        MF1 <- tryCatch(
          {
            # Denominator model (YEAR + treatment(reporting) history variables + time-varying covariates)
            denom.fit <-  glm(CATreview ~ lngdpl1 + lnpop 
                              +nhril1 +polityl1+ tortl1 + HRTpercentl1 + CATart22l1
                              +   CATrevprevyr + CATyrssincereview + CATnumprevrevs+nevdem +
                                CATregiondensityl1+regionall1+transdeml1+ transition2l1+ yr,
                              family="binomial", data=boot.data)
            denom.p <- predict(denom.fit, type = "response")
            
            # Numerator model (YEAR + treatment(reporting) history variables)
            numer.fit <-  glm(CATreview ~ 
                                CATrevprevyr + CATyrssincereview + CATnumprevrevs + 
                                nevdem + yr, 
                              family="binomial", data=boot.data)
            numer.p <- predict(numer.fit, type = "response")
            if (denom.fit$converged == 0 || numer.fit$converged == 0){
              stop("Probability models didn't converge")
              
            }else{
              ## Calculate weights
              boot.data$treat.prob <- ifelse(boot.data$CATreview == 1, denom.p, denom.p)
              boot.data$d.pr.tr <- ifelse(boot.data$CATreview == 1, denom.p, 1-denom.p)
              boot.data$n.pr.tr <- ifelse(boot.data$CATreview == 1, numer.p, 1-numer.p)
              boot.data$w <- boot.data$n.pr.tr/boot.data$d.pr.tr
              
              ## Bias-adjusted outcome - alpha*(1-1)*Pr(1) + alpha*(1-0)*Pr(0) for units receiving treatment and alpha*(0-1)*Pr(1) + alpha(0-0)*Pr(0) for units receiving control
              boot.data$time.bias <- ifelse(boot.data$CATreview == 1, alpha*(1-boot.data$treat.prob), -alpha*(boot.data$treat.prob))
              
              ## Observation weight is product over time
              boot.data.dt <- data.table(boot.data)
              # Product of weights
              boot.data.dt[ intersect(order(year), which(!is.na(w))), swrev := cumprod(w), by=uniqueID]
              # Bias at time period
              boot.data.dt[ intersect(order(year), which(!is.na(time.bias))), s.time.bias := cumsum(time.bias), by=uniqueID]
              boot.data <- as.data.frame(boot.data.dt)
              
              
              
              ###################################################
              ##### Sensitivity analysis transformations
              ###################################################
              if (alpha !=0){
                boot.data$outcomeAdj <- boot.data$outcome - boot.data$s.time.bias
              }else{
                boot.data$outcomeAdj <-factor(boot.data$outcome,ordered=TRUE)
              }
              
              ########
              polr(outcomeAdj ~ CATreview + CATrevprevyr+ CATyrssincereview+ 
                     CATnumprevrevs+ CATlintt ,method='probit',  weights = boot.data$swrev, 
                   data = boot.data)
            }
          }, error = function(cond) {
            
            return("fail")
          }
        )
        ## If no error was caught
        if (class(MF1) == "polr"){
          ### Save the coefficients
          boot_coefs[iter,] <-c(MF1$coefficients,MF1$zeta)
          if (iter == 1){
            colnames(boot_coefs) <-c(names(MF1$coefficients), names(MF1$zeta))
          }
          
          ## Increment iterations by 1
          iter <- iter + 1
        }else{
          print("Error detected! Re-sampling")
        }
      }
      
      #### Review model  
    }else if (model == 2){
      ## Calculate number of country clusters
      n_country <- length(unique(catuse$cowcode))
      
      ## Number of coefficients
      n_coef <- 5
      
      ## Placeholder for holding bootstrap coefficients
      boot_coefs <- as.data.frame(matrix(nrow=iterations, ncol=(n_coef+2)))
      
      # Start bootstrap loop - since we're concerned about sampling datasets that have non-convergence
      # we make this a while loop and use tryCatch when dealing with non-feasible datasets.
      
      # Initialize iteration number
      iter <- 1
      while(iter <= iterations){
        print(iter)
        # If it's not the first run, bootstrap
        if (iter != 1){
          # Resample dataset by country clusters
          boot.data <- data.frame()
          for (kit in 1:n_country){
            cnum <- sample(unique(catuse$cowcode), size=1) # Sample a country
            sampled.country <- catuse[catuse$cowcode == cnum,]
            sampled.country$uniqueID <- kit
            boot.data <- rbind(boot.data, sampled.country) # Append to the bootstrap dataset
          }
          # If it's the first run, use the raw data
        }else{
          boot.data <- catuse
          boot.data$uniqueID <- boot.data$cowcode
        }
        
        ### Get Coefficients - catch any errors
        MF1 <- tryCatch(
          {
            # Denominator model (YEAR + treatment(review) history variables + time-varying covariates)
            denom.fit <-  glm(CATreview ~ lngdpl1 + lnpop 
                              +nhril1 +polityl1+ tortl1 + HRTpercentl1 + CATart22l1
                              +   CATrevprevyr + CATyrssincereview + CATnumprevrevs+nevdem +
                                CATregiondensityl1+regionall1+transdeml1+transition2l1+ transparencyindex_l1+
                                yr, 
                              family="binomial", data=boot.data)
            denom.p <- predict(denom.fit, type = "response")
            
            # Numerator model (YEAR + treatment(review) history variables)
            numer.fit <-  glm(CATreview ~ 
                                CATrevprevyr + CATyrssincereview + CATnumprevrevs + nevdem +
                                yr, 
                              family="binomial", data=boot.data)
            numer.p <- predict(numer.fit, type = "response")
            
            
            if (denom.fit$converged == 0 || numer.fit$converged == 0){
              stop("Probability models didn't converge")
            }else{
              ## Calculate weights
              boot.data$treat.prob <- ifelse(boot.data$CATreview == 1, denom.p, denom.p)
              boot.data$d.pr.tr <- ifelse(boot.data$CATreview == 1, denom.p, 1-denom.p)
              boot.data$n.pr.tr <- ifelse(boot.data$CATreview == 1, numer.p, 1-numer.p)
              boot.data$w <- boot.data$n.pr.tr/boot.data$d.pr.tr
              
              ## Bias-adjusted outcome - alpha*(1-1)*Pr(1) + alpha*(1-0)*Pr(0) for units receiving treatment and alpha*(0-1)*Pr(1) + alpha(0-0)*Pr(0) for units receiving control
              boot.data$time.bias <- ifelse(boot.data$CATreview == 1, alpha*(1-boot.data$treat.prob), -alpha*(boot.data$treat.prob))
              
              ## Observation weight is product over time
              boot.data.dt <- data.table(boot.data)
              # Product of weights
              boot.data.dt[ intersect(order(year), which(!is.na(w))), swrev := cumprod(w), by=uniqueID]
              # Bias at time period
              boot.data.dt[ intersect(order(year), which(!is.na(time.bias))), s.time.bias := cumsum(time.bias), by=uniqueID]
              boot.data <- as.data.frame(boot.data.dt)
              
              
              ###################################################
              ##### Sensitivity analysis transformations
              ###################################################
              if (alpha !=0){
                boot.data$outcomeAdj <- boot.data$outcome - boot.data$s.time.bias
              }else{
                boot.data$outcomeAdj <-factor(boot.data$outcome,ordered=TRUE)
              }
              
              ########
              ## Estimate MSM
              polr(outcomeAdj ~ CATreview + CATrevprevyr+ CATyrssincereview+ 
                     CATnumprevrevs+ CATlintt ,method='probit',  weights = boot.data$swrev, 
                   data = boot.data)
            }
          }, error = function(cond) {
            
            return("fail")
          }
        )
        ## If no error was caught
        if (class(MF1) == "polr"){
          ### Save the coefficients
          boot_coefs[iter,] <-c(MF1$coefficients,MF1$zeta)
          if (iter == 1){
            colnames(boot_coefs) <-c(names(MF1$coefficients), names(MF1$zeta))
          }
          
          ## Increment iterations by 1
          iter <- iter + 1
        }else{
          print("Error detected! Re-sampling")
        }
        
      }
      ### Model 3: Interaction model
    }else if (model == 3){
      ## Calculate number of country clusters
      n_country <- length(unique(catuse$cowcode))
      
      ## Number of coefficients
      n_coef <- 7
      
      ## Placeholder for holding bootstrap coefficients
      boot_coefs <- as.data.frame(matrix(nrow=iterations, ncol=(n_coef+2)))
      
      # Start bootstrap loop - since we're concerned about sampling datasets that have non-convergence
      # we make this a while loop and use tryCatch when dealing with non-feasible datasets.
      
      # Initialize iteration number
      iter <- 1
      while(iter <= iterations){
        print(iter)
        # If it's not the first run, bootstrap
        if (iter != 1){
          # Resample dataset by country clusters
          boot.data <- data.frame()
          for (kit in 1:n_country){
            cnum <- sample(unique(catuse$cowcode), size=1) # Sample a country
            sampled.country <- catuse[catuse$cowcode == cnum,]
            sampled.country$uniqueID <- kit
            boot.data <- rbind(boot.data, sampled.country) # Append to the bootstrap dataset
          }
          # If it's the first run, use the raw data
        }else{
          boot.data <- catuse
          boot.data$uniqueID <- boot.data$cowcode
        }
        
        ### Get Coefficients - catch any errors
        MF1 <- tryCatch(
          {
            # Denominator model (YEAR + treatment(review) history variables + time-varying covariates)
            denom.fit <-  glm(CATreview ~ lngdpl1 + lnpop 
                              +nhril1 +polityl1+ tortl1 +farissl2+ HRTpercentl1 + CATart22l1
                              +   CATrevprevyr + CATyrssincereview + CATnumprevrevs+nevdem +
                                CATregiondensityl1+regionall1+transdeml1+transition2l1+ 
                                transparencyindex_l1 +CATlintt+farissl2*CATlintt, 
                              family="binomial", data=boot.data)
            denom.p <- predict(denom.fit, type = "response")
            
            # Numerator model (YEAR + treatment(review) history variables)
            numer.fit <-  glm(CATreview ~ 
                                CATrevprevyr + CATyrssincereview + CATnumprevrevs + nevdem+ CATlintt, 
                              family="binomial", data=boot.data)
            numer.p <- predict(numer.fit, type = "response")
            if (denom.fit$converged == 0 || numer.fit$converged == 0){
              stop("Probability models didn't converge")
            }else{
              ## Calculate weights
              boot.data$treat.prob <- ifelse(boot.data$CATreview == 1, denom.p, denom.p)
              boot.data$d.pr.tr <- ifelse(boot.data$CATreview == 1, denom.p, 1-denom.p)
              boot.data$n.pr.tr <- ifelse(boot.data$CATreview == 1, numer.p, 1-numer.p)
              boot.data$w <- boot.data$n.pr.tr/boot.data$d.pr.tr
              
              ## Bias-adjusted outcome - alpha*(1-1)*Pr(1) + alpha*(1-0)*Pr(0) for units receiving treatment and alpha*(0-1)*Pr(1) + alpha(0-0)*Pr(0) for units receiving control
              boot.data$time.bias <- ifelse(boot.data$CATreview == 1, alpha*(1-boot.data$treat.prob), -alpha*(boot.data$treat.prob))
              
              ## Observation weight is product over time
              boot.data.dt <- data.table(boot.data)
              # Product of weights
              boot.data.dt[ intersect(order(year), which(!is.na(w))), swrev := cumprod(w), by=uniqueID]
              # Bias at time period
              boot.data.dt[ intersect(order(year), which(!is.na(time.bias))), s.time.bias := cumsum(time.bias), by=uniqueID]
              boot.data <- as.data.frame(boot.data.dt)
              
              
              ###################################################
              ##### Sensitivity analysis transformations
              ###################################################
              if (alpha !=0){
                boot.data$outcomeAdj <- boot.data$outcome - boot.data$s.time.bias
              }else{
                boot.data$outcomeAdj <-factor(boot.data$outcome,ordered=TRUE)
              }
              
              
              ########
              ## Estimate MSM
              polr(outcomeAdj ~ CATreview + CATrevprevyr+ CATyrssincereview+ 
                     CATnumprevrevs+ farissf1+CATlintt+farissf1*CATlintt ,method='probit',  weights = boot.data$swrev, 
                   data = boot.data)
            }
          }, error = function(cond) {
            
            return("fail")
          }
        )
        
        ## If no error was caught
        if (class(MF1) == "polr"){
          ### Save the coefficients
          boot_coefs[iter,] <-c(MF1$coefficients,MF1$zeta)
          if (iter == 1){
            colnames(boot_coefs) <-c(names(MF1$coefficients), names(MF1$zeta))
          }
          
          ## Increment iterations by 1
          iter <- iter + 1
        }else{
          print("Error detected! Re-sampling")
        }
        
      }
      
    }
    
    boot_coefs$model <- model
    boot_coefs$alpha <- alpha
    boot_coefs$iterations <- iterations
    boot_coefs$outcomelead <- lead
    ## Save results
    filen <- paste("CAT_CIRI__model_",model,"_alpha_", alpha, "_iterations_", iterations, "_outcomelead_",lead, ".csv", sep="")
    write.csv(boot_coefs, filen)
    sink(NULL)
  }
}

## Get the results directory

sink("CAT CIRI MSM Analysis.txt",append=TRUE)
## Load Models - base results

model1<-read.csv("CAT_CIRI__model_1_alpha_0_iterations_20000_outcomelead_2.csv", header=TRUE)
m1_coeff<-rbind(model1$CATreview,model1$CATrevprevyr,model1$CATyrssincereview,model1$CATnumprevrevs,model1$CATlintt)

####Code to obtain the mean coefficients and their 95% CI
m1_95ci_upper <- apply(model1, 2, function(x) quantile(x, .975))
m1_95ci_lower <- apply(model1, 2, function(x) quantile(x, .025))
m1_mean<- apply(model1, 2, function(x) quantile(x, .50))
out1mean<-(m1_mean)
out1low<-(m1_95ci_lower)
out1upper<-(m1_95ci_upper)
g1<-rbind(out1mean,out1low,out1upper)
ci1<-t(g1)

stargazer(ci1)

xc.nrrev<-cbind(0, 0, 2.5, c(0:5), 13.5)
head(xc.nrrev)

ystar.nrrev<-(xc.nrrev)%*%(m1_coeff)
ystar.nrrev<-t(ystar.nrrev)
dim(ystar.nrrev)
m1_zeta<-cbind(model1$X0.1, model1$X1.2)
m1_zeta<-as.data.frame(m1_zeta)
###
pr1<- pnorm(m1_zeta$V1,mean=ystar.nrrev,sd=1)
pr2<- pnorm(m1_zeta$V2,mean=ystar.nrrev,sd=1)-pnorm(m1_zeta$V1,mean=ystar.nrrev,sd=1)
pr3<-1-pnorm(m1_zeta$V2,mean=ystar.nrrev,sd=1)

pr1mean<-apply(pr1,2,median)
pr2mean<-apply(pr2,2,median)
pr3mean<-apply(pr3,2,median)




postscript("FigureA5.eps", width=4.5, height=4.5, horizontal = FALSE, 
           onefile = FALSE)
par(cex.axis=1.5, cex.lab=1.2, cex.main=1.7, cex.sub=1.7)
plot(c(0:5),pr1mean,ylim=c(0,1),xlim=c(0,5),lwd=2,type="l",lty=3,
     xlab="Number of Previous Reviews ", ylab="Probability of observing a given level of Torture")

lines(c(0:5),pr2mean,ylim=c(0,1),lwd=2,lty=2,type="l")
lines(c(0:5),pr3mean,ylim=c(0,1),lwd=2,lty=4,type="l")
legend(1,1, c("No torture","Some torture","Widespread torture"),
       lty=c(3,0,0),lwd=2, bty="n",cex=.9)
legend(1,1, c("No torture","Some torture","Widespread torture"),
       lty=c(0,2,0),lwd=2, bty="n",cex=.9)
legend(1,1, c("No torture","Some torture","Widespread torture"),
       lty=c(0,0,4),lwd=2, bty="n",cex=.9)

p1.low<- apply(pr1,2,quantile,prob=.025)
p1.high<- apply(pr1,2,quantile,prob=.975)
segments(x0=seq(0,5,1), 
         y0=p1.low,  
         y1=p1.high,lwd=1) 
p2.low<- apply(pr2,2,quantile,prob=.025)
p2.high<- apply(pr2,2,quantile,prob=.975)
segments(x0=seq(0,5,1), 
         y0=p2.low,  
         y1=p2.high,lwd=1) 
p3.low<- apply(pr3,2,quantile,prob=.025)
p3.high<- apply(pr3,2,quantile,prob=.975)
segments(x0=seq(0,5,1), 
         y0=p3.low,  
         y1=p3.high,lwd=1)
dev.off()


########Figure A6
model2<-read.csv("CAT_CIRI__model_2_alpha_0_iterations_20000_outcomelead_2.csv")
m2_coeff<-rbind(model2$CATreview,model2$CATrevprevyr,model2$CATyrssincereview,model2$CATnumprevrevs,model2$CATlintt)
m2_95ci_upper <- apply(model2, 2, function(x) quantile(x, .975))
m2_95ci_lower <- apply(model2, 2, function(x) quantile(x, .025))
m2_mean<- apply(model2, 2, function(x) quantile(x, .50))

out2mean<-(m2_mean)
out2low<-(m2_95ci_lower)
out2upper<-(m2_95ci_upper)
g2<-rbind(out2mean,out2low,out2upper)
ci2<-t(g2)
stargazer(ci2)

xc.nrrev<-cbind(0, 0, 2.63, c(0:5), 13.2)
head(xc.nrrev)

ystar.nrrev<-(xc.nrrev)%*%(m2_coeff)
ystar.nrrev<-t(ystar.nrrev)
dim(ystar.nrrev)
m2_zeta<-cbind(model2$X0.1, model2$X1.2)
m2_zeta<-as.data.frame(m2_zeta)
###
pr1<- pnorm(m2_zeta$V1,mean=ystar.nrrev,sd=1)
pr2<- pnorm(m2_zeta$V2,mean=ystar.nrrev,sd=1)-pnorm(m2_zeta$V1,mean=ystar.nrrev,sd=1)
pr3<-1-pnorm(m2_zeta$V2,mean=ystar.nrrev,sd=1)
pr1mean<-apply(pr1,2,median)
pr2mean<-apply(pr2,2,median)
pr3mean<-apply(pr3,2,median)



postscript("FigureA6.eps", width=4.5, height=4.5, horizontal = FALSE, 
           onefile = FALSE)
par(cex.axis=1.5, cex.lab=1.2, cex.main=1.7, cex.sub=1.7)
plot(c(0:5),pr1mean,ylim=c(0,1),xlim=c(0,5),lwd=2,type="l",lty=3,
     xlab="Number of Previous Reviews ", ylab="Probability of observing a given level of Torture")

lines(c(0:5),pr2mean,ylim=c(0,1),lwd=2,lty=2,type="l")
lines(c(0:5),pr3mean,ylim=c(0,1),lwd=2,lty=4,type="l")
legend(1,1, c("No torture","Some torture","Widespread torture"),
       lty=c(3,0,0),lwd=2, bty="n",cex=.9)
legend(1,1, c("No torture","Some torture","Widespread torture"),
       lty=c(0,2,0),lwd=2, bty="n",cex=.9)
legend(1,1, c("No torture","Some torture","Widespread torture"),
       lty=c(0,0,4),lwd=2, bty="n",cex=.9)

p1.low<- apply(pr1,2,quantile,prob=.025)
p1.high<- apply(pr1,2,quantile,prob=.975)
segments(x0=seq(0,5,1), 
         y0=p1.low,  
         y1=p1.high,lwd=1) 
p2.low<- apply(pr2,2,quantile,prob=.025)
p2.high<- apply(pr2,2,quantile,prob=.975)
segments(x0=seq(0,5,1), 
         y0=p2.low,  
         y1=p2.high,lwd=1) 
p3.low<- apply(pr3,2,quantile,prob=.025)
p3.high<- apply(pr3,2,quantile,prob=.975)
segments(x0=seq(0,5,1), 
         y0=p3.low,  
         y1=p3.high,lwd=1)
dev.off()


sink(NULL)
